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Thermoelectric  devices  are  semiconductor  devices  which  are  capable  of  either  generating  a  voltage  when 
placed  in  between  a  temperature  gradient,  exploiting  the  Seebeck  effect,  or  producing  a  temperature 
gradient  when  powered  by  electricity,  exploiting  the  Peltier  effect.  The  devices  are  usually  employed  in 
environments  with  time-varying  temperature  differences  and  input/output  powers.  Therefore  it 
becomes  important  to  understand  the  behaviour  of  thermoelectric  devices  during  thermal  and  electrical 
transients  in  order  to  properly  simulate  and  design  complex  thermoelectric  systems  which  also  include 
power  electronics  and  control  systems. 

The  purpose  of  this  paper  is  to  provide  the  transient  solution  to  the  one-dimensional  heat  conduction 
equation  with  internal  heat  generation  that  describes  the  transfer  and  generation  of  heat  throughout 
a  thermoelectric  device.  The  solution  proposed  can  be  included  in  a  model  in  which  the  Peltier  effect,  the 
thermal  masses  and  the  electrical  behaviour  of  the  system  are  considered  too;  this  would  be  of  great 
benefit  because  it  would  allow  accurate  simulations  of  thermoelectric  systems. 

While  the  previous  literature  does  not  focus  on  the  study  of  thermal  transients  in  thermoelectric 
applications  and  usually  considers  constant  the  temperatures  at  the  hot  and  cold  sides,  this  paper 
proposes  a  dynamic  exchange  of  heat  through  the  hot  and  cold  side,  both  in  steady-state  and  transients. 
This  paper  also  presents  an  analytical  solution  which  is  then  computed  by  Matlab  to  simulate  a  physical 
experiment.  Simulation  results  show  excellent  correlation  with  experimentally  determined  values,  thus 
validating  the  solution. 

©  2011  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

Current  literature  [1,2]  considers  thermoelectric  (TE)  systems 
under  steady-state  and  when  in  equilibrium.  Few  papers  [3-7] 
consider  the  simulation  of  thermal  transients  found  in  practical 
TE  systems,  or  effects  due  to  variation  in  the  electrical  system 
parameters.  Actually  most  thermoelectric  applications  are  subject 
to  electro-thermal  transients.  In  cooling  applications  for  electronic 
devices  [8]  heat  is  generated  depending  on  use,  therefore 
a  controller  is  usually  employed  and  its  design  needs  to  take 
transients  into  account.  In  thermoelectric  power  generation 
temperatures  are  often  changeable,  especially  in  applications  to  the 
automotive  field  [9,10]  or  to  stoves  [11,12],  where  start-up  and  shut¬ 
down  considerations  are  a  major  concern.  These  problems  are  yet 
not  fully  assessed  in  literature  and  industry. 

The  following  part  of  this  section  focuses  on  TE  generators 
(TEGs),  but  similar  considerations  can  be  applied  to  Peltier  devices. 
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The  operation  of  thermoelectric  devices  depends  on  five  main 
effects:  (1)  thermal  conduction;  (2)  internal  Joule  heating;  (3) 
Seebeck  effect;  (4)  Peltier  effect  and  (5)  Thomson  phenomenon, 
which  is  here  neglected  because  of  small  effect  [1,2]. 

The  thermal  power  input  to  the  hot  junction  is  given  by 

MAT  1  2 

Qh  =  — j - f  aTHI  -  ^R[ntr  (1 ) 

where  k  is  the  overall  conduction  coefficient,  A  is  the  area  and  L  the 
thickness  of  the  TEG,  AT  is  the  temperature  gradient,  a  is  the  See¬ 
beck  coefficient,  TH  is  the  temperature  at  the  hot  side,  I  is  the 
current  produced  and  R[nt  is  the  overall  internal  resistance  of  the 
device. 

As  explained  in  a  paper  by  Ben-Yaakov  and  Lineykin  [3],  Eq.  (1 )  is 
derived  from  the  steady-state  solution  of  the  one-dimensional  heat 
conduction  equation  for  solids  with  internal  energy  generation, 
which  is  written  as: 
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where  g  is  the  rate  of  heat  internally  generated  per  unit  volume  [W/ 
m3]  while  T  and  k  have  the  same  meaning  as  in  Eq.  (1). 

Assuming  constant  temperatures  at  the  hot  and  cold  side  as 
boundary  conditions  (T(t, 0)  =  TH;  T(t,L )  =  Tc),  Eq.  (2)  can  be  solved  as: 


T  W.T„-!*2 


Tc-Th 


—L2 
2k  , 


(3) 


Combining  this  result  to  the  equation  of  heat  transfer 
q  =  (-/<A)ar/ax,  the  expression  for  the  heat  flow  through  the  hot 
side  becomes 


Qh 


/<AAT 

L 


~2Rintl2 


(4) 


Taking  into  account  the  Peltier  heating  or  cooling  effect  at  the 
junctions  this  equation  becomes  Eq.  (1). 

The  solution  for  the  cold  side  is  similar;  the  two  results  show 
that  the  Joule  heating  in  the  bulk  material  is  equally  divided 
between  the  hot  and  cold  sides.  This  steady-state  characterisation 
of  TE  devices  accurately  describes  the  balance  of  powers  in  a  system 
under  steady-state.  However  it  is  not  usable  in  a  time-varying 
environment  because  it  does  not  take  time  into  consideration  and 
because  the  boundary  conditions  used  to  solve  Eq.  (2)  assume  that 
the  temperatures  at  the  hot  and  cold  sides  are  constant. 

Test  results  from  a  paper  by  Chen  L.  et  al.  [4],  which  investigates 
the  electrical  response  of  TE  devices  to  load  transients,  show  that 
a  common  TE  module  has  a  very  fast  dynamic  response,  in  the  order 
of  nano-seconds.  This  means  that  the  direct  energy  conversion 
taking  place  in  the  practical  TE  systems  can  be  considered  instan¬ 
taneous.  Hence  it  does  not  influence  the  thermal  dynamics,  which 
are  controlled  only  by  the  thermal  elements,  i.e.  thermal  capacities 
of  the  system. 

The  thermal  time  constants  of  TE  devices  are  of  higher  orders  of 
magnitude  if  compared  to  the  electrical  ones.  This  means  that  there 
are  “spikes”  in  power  and  efficiency  before  the  thermal  time 
constants  bring  the  device  to  steady-state  operation.  When  the 
operating  conditions  change  simultaneously  with  time,  the  effects 
of  these  transients  are  so  important  that  both  the  heat  sources  and 
the  electronics  is  affected.  Therefore  their  impact  on  the  entire 
system,  besides  the  steady-state  behaviour,  is  important.  Depend¬ 
ing  on  the  power  input  to  the  hot  side,  the  number  of  modules  and 
the  power  output  from  the  cold  side,  these  thermal  transients  can 
last  for  quite  long  periods  of  time,  in  the  range  of  tens  of  seconds.  To 
carefully  design  the  whole  system  in  cases  of  frequent  changes  of 
operation  conditions  and  during  start-up  and  shut-down,  it 
becomes  necessary  to  model  these  thermal  transients  and  dynamic 
characteristics. 

In  both  papers  by  Chen  M.  et  al.  [5],  as  well  as  in  [3],  the  transient 
term  is  included  in  the  (SPICE)  model  as  a  parallel  electrical 
capacitor,  converting  the  thermal  model  into  an  equivalent  elec¬ 
trical  circuit,  to  take  into  account  an  additional  time-related  term 
due  to  the  change  in  stored  heat  energy.  Therefore,  considering  an 
infinitesimal  element  dx  inside  the  TE  module,  as  in  Fig.  1,  the  heat 
flow  increase  from  the  inlet  to  the  outlet  of  the  element  is 

d q  =  q(x  +  dx)  -  q(x)  =  RintI2  +  mCm ^  (5) 

where  m  and  Cm  are  the  mass  and  the  specific  heat  of  the  element. 
0  is  the  thermal  resistance  of  dx. 

A  commonly  used  approach  is  to  divide  the  module  into  a  grid  of 
small  elements  like  dx,  each  one  taking  into  account  the  heat 
conduction,  the  Joule  heating,  the  Peltier  heating/cooling  and  the 
thermal  mass;  hence  an  accurate  equivalent  circuit  is  composed  by 
a  series  of  RC  cells.  Additionally,  other  parts  of  the  TE  system,  like 
the  ceramic  plates  and  the  heat  sinks,  need  to  be  included  as 


dx 


Fig.  1.  Lumped  equivalent  circuit  of  the  heat  transfer  in  an  infinitesimal  volume  dx 
inside  the  TE  module  [5]. 

thermal  impedances.  The  accuracy  of  the  simulation  is  related  to 
the  number  of  cells  used,  which  needs  to  be  tuned  depending  on  if 
the  thermal  masses  of  these  external  elements  are  comparable, 
bigger  or  smaller  than  the  one  of  the  TE  module.  These  kinds  of 
models  available  for  circuit  simulators  like  SPICE  and  ANSYS  do  not 
offer  a  theoretical  solution  to  the  problem  and  are  sometimes 
difficult  to  use,  due  to  the  division  of  the  module  in  multiple 
elements,  whose  parameters  are  difficult  to  obtain  from  the  data 
sheets  of  the  manufacturers. 

A  more  appropriate  and  more  reliable  approach  would  be  to 
study  these  devices  through  the  physical  equations  which  describe 
their  behaviour;  among  these,  the  most  important  is  the  heat 
equation.  Alata,  Naji  and  Al-Nimr  [6,7]  have  already  explored  this 
possibility,  but  they  used  fixed  temperatures  as  boundary  condi¬ 
tions  at  the  two  sides  of  the  TE  module,  assuming  that  those 
temperatures  are  not  varying  even  if  there  is  an  exchange  of  heat 
through  the  sides,  i.e.  supposing  thermal  isolation. 

The  following  section  describes  a  mathematical  solution  of  the 
heat  conduction  equation  for  TE  devices,  with  internal  Joule  heat 
generation,  and  dynamic  exchanges  of  heat  through  the  hot  and 
cold  sides. 


2.  Theory  and  calculation 


In  TE  devices  the  Peltier  heating  or  cooling  acts  only  at  the 
junctions  where  there  is  the  connection  with  metal,  therefore, 
assuming  there  are  no  temperature  gradients  and  losses  in  the  y 
and  z  directions,  the  heat  propagation  and  generation  inside  the 
module  can  be  modelled  through  the  one-dimensional,  unsteady 
heat  conduction  equation  with  constant  thermal  conductivity  and 
internal  heat  generation,  expressed  as 


92T  g  _  1  9T 
9x2  k  £  dt 


(6) 


where  T  is  the  temperature  in  Kelvin  degrees  function  of  both  the 
space  x  and  the  time  t,  g  is  the  rate  of  heat  energy  generation  per 
unit  volume  [W/m3],  k  is  the  thermal  conductivity  coefficient  [W/ 
m  K],  and  e  is  the  thermal  diffusivity  coefficient  defined  as 


_  heat  conducted  _  k 
~  heat  stored  —  pCm 

where  pCm  is  the  heat  capacity  per  unit  volume  [J/m3  K]. 

Our  objective  is  to  study  TE  devices  under  conditions  where 
both  the  time  and  temperatures  are  varying,  hence  the  tempera¬ 
tures  at  the  hot  and  cold  sides  are  not  considered  constant.  They  are 
on  the  contrary  considered  connected  to  ’’far-away”  constant 
temperatures,  thus  allowing  their  temperatures  to  vary  depending 
on  the  rate  of  heat  flow  through  them  and  the  generation  of  heat 
inside  the  TE  device. 
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2.1.  Description  of  the  problem 

To  simplify  the  symmetry  of  the  problem,  consider  a  TE  module 
of  thickness  2 L  and  area  A ,  positioned  on  the  x-axis  in  such  a  way 
that  the  hot  side  is  at  x  =  -L  and  the  cold  side  at  x  =  L.  Hence, 
simplifying  the  notation,  consider  the  following  equation: 

Tm+i  =  -Tt  (8) 

l<  e 

on  the  interval  [— L,L]. 

In  order  to  be  able  to  solve  this  differential  differential  equation 
it  is  necessary  to  have  two  boundary  conditions  and  an  initial  state. 
The  latter  is  assumed  as  a  linear  distribution  of  temperature 
throughout  the  whole  module,  starting  from  the  initial  hot 
temperature  Tni  atx  =  -L,  to  the  initial  cold  temperature  To  atx  =  L. 
Hence  the  initial  condition  is 


%  +  y  =  Ox 

and  the  boundary  conditions  in  Eq.  (10)  become 

(12) 

0t  =  fa(6- 1),  at£  =  -l 

(13a) 

0f  =  -0c(0+ 1).  at?  =  l 

(13b) 

Finally,  the  initial  state  in  Eq.  (9)  becomes 

=  (Tci  ~  Th£  ~  (JH oq  +  TCoo  -  THi  -  Ta)) 

TH  00  —  TCo o 

The  solution  of  this  problem  will  be  the  sum  of  a  steady-state 
solution  and  a  transient  solution,  as  described  in  the  next  two 
sub-sections.  The  final  solution  for  6  can  then  be  written  as 


T(x,0)  =  T0(x)  -  (9) 

The  boundary  conditions  are  those  appropriate  to  the  Newton’s 
Law  of  Cooling,  stating  that  the  heat  flux  (per  unit  area)  out  of 
a  boundary  with  normal  n  is  q  =  -kn  -VT,  so  that  the  boundary 
conditions  for  this  problem  are  written  as 


0(?,T)  =  0SS(f)*0t(f,T)  (15) 

where  dss  is  the  steady-state  profile  and  dt  is  the  transient 
evolution. 

2.3.  Steady-state  solution 


Tx  =  @h(T  —  THoo),  atx  =  -L  (10a) 

Tx  =  -VC(J-TC  „),  atx  =  L  (10b) 

where  Th™  and  Tc «>  are  the  temperatures  “far-away”  from  the  hot 
and  cold  sides  of  the  TEG.  Fig.  2  illustrates  the  physical  system.  (3H 
and  (3C  are  constants. 


2.2.  Scaling 


It  is  convenient  to  appropriately  scale  the  problem  in  order  to 
identify  the  important  parameter  combinations: 


L 2 

x  =  L?;  t  =  — T 
£ 

T  =  AT0  +  —  +  — 


Ph:c 


h±y  =  #L- 

L  1  J  /<AT 


(ID 


where  AT  =  {TH «>  -  TCoo)/2.  The  parameter  y  represents  how  much 
heat  is  internally  generated  by  Joule  heating  compared  to  the  heat 
flux  through  the  device  for  conduction,  and  this  is  due  to  the 
temperature  gradient  across  the  device. 

After  this  scaling  the  problem  in  Eq.  (8)  becomes 


We  wish  to  find  the  steady-state  part  of  the  solution  in  Eq.  (15). 
If  t  — >  oo  then  a/a t  ->  0,  that  is  Ox  =  0.  Hence  the  Poisson’s  equation 
hi  +  T  =  0,  which  leads  to 

tfss(?)  +  (16) 

where  A  and  B  must  be  chosen  to  satisfy  the  boundary  conditions, 


namely: 

B  +  y  =  0h(A-B-l- l) 


B-y  =  — 0c(a  +  b-|  +  i) 

leading  to 

T  ^2  +  —fa  +  —@c  +  PcPh^j 
Pc  +  Ph  +  cPh 

p  =  T  (Ph-Pc)-2PcPh 
Pc  +  Ph  +  2PcPjj 


(17a) 

(17b) 

Ph-Pc 

-  (18a) 

(18b) 


2.4.  Transient  solution 

The  transient  part  of  Eq.  (15),  6t  is 


%  *  h 

(19) 

with  boundary  conditions 

-  ph0  =  0  at  f  =  -1 

(20a) 

h  +  Pc0  =  0  at  ?  =  1 

(20b) 

Eq.  (19)  is  the  heat  equation  with  homogeneous  boundary 
conditions  (Eq.  (20));  as  such  we  search  for  separable  solutions  and 
using  the  transform  methods  (Fourier  transform  on  finite  domain 
gives  Fourier  series)  we  arrive  at  the  following  expression  for  dt: 
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(Akcosk%  +  BksinkZ)e  k2z  (21) 

allowable  k 


where  the  coefficients  Ak  and  Bk  and  the  allowable  wave-numbers  k 
are  determined  by  the  boundary  conditions  and  the  initial  condi¬ 
tion.  First  we  consider  the  allowable  k.  The  boundary  conditions  of 
Eq.  (20)  give 


i 

J  V'k(?)V'k(?)d?  =cos2k[/c2  -&A 
-1 


(29) 


+  sin2k 


5  +  TF  +  W  +  W 


Suppose  the  initial  condition  is  0q(?).  then  at  t  =  0,  we  have  that 

0t(?,O)  =  0o(f)-Mf).  Defining 


kA^sink  +  kBkcosk  -  fih(Akcosk  -  Bksink)  =  0  (22a) 

-kAksink  +  kBkcosk  +  Pc(Akcosk  +  6/<sink)  =  0  (22b) 

which  can  be  written  in  matrix  form  as 


i6ccos k  -  ksink 

i6csink  +  kcosk 

O' 

-fihcosk  +  ksink 

ifyjSink  +  kcosk 

lBk\ 

0 

whose  determinant  is 

|det|  =  (ficfih  ~  k2)  sin2/<  +  k((ic  +  /?h)cos2 k  (23) 

which  has  non-trivial  solutions  (|det|=0)  for  values  of  k  such  that 


l 

/  ^(f)(«o(f)-»ss(f))df 

4.  =  - - T -  (30) 

-i 

Then  the  solution  for  6  is 

KM  =  OssM  +  n'EKMe-*  (3D 

The  final  solution  in  Eq.  (31 )  obviously  needs  to  be  scaled  back  to 
temperature  using  the  reverse  scaling  of  Eq.  (11). 


tan2k  = 


+  Ph) 

k2  -  Pch 


(24) 


To  find  a  corresponding  eigenfunction  let’s  consider  Eq.  (22a)  in 
the  form 


Afc(ksink  -  (3hcos k)  +  6fc(kcosk  +  (3hsink)  =  0  (25) 

A  corresponding  eigenfunction  is 


2.6.  The  boundary  conditions 

The  boundary  conditions  used  to  solve  the  heat  equation  come 
from  the  Newton’s  Law  of  Cooling  and  they  are  written  in  the  form 
of  Eq.  (10).  It  is  now  important  to  define  the  coefficient  /3[m— 1  ]  in 
physical  terms. 

Considering  a  body  of  volume  B,  with  boundaries  66,  the  total 
energy  in  the  body  is 


<Mf)  =  -0csin fe(?  -  1)  +  fccosfc(?  -  1)  (26) 

Enumerating  the  solutions  of  Eq.  (24)  as  kn,  with 
0  </<!<•••<  kn  <  •  •  •  for  ne  N,  it  can  be  noted  that  tan2kn  =  0  for  big 
values  of  /<n.  This  happens  when  2 kn  =  mru,  with  m  integer,  that  is 
when  kn  =  nuz/2. 

The  solution  for  0t  can  then  be  written  as 

KiM  =  'Z'K'hc.e-*  (27) 

n 

where  the  coefficients  A'n  are  determined  by  the  initial  conditions 
0o(?)  of  Eq.  (14),  as  explained  in  the  following  section. 


2.5.  Transforming  the  initial  condition 

The  operator  defined  by  Lk6  =  6%  +  k2d  =  0  with  the  boundary 
conditions  of  Eq.  (20)  is  self-adjoint. 

Some  notation,  let 


E  =  J  pCmTdV  (32) 

B 

where  pCm  represents  the  heat  storage  capability  of  the  body  per 
unit  volume  [J/m3  I<]. 

From  the  heat  flux  q  =  -]<■  VT and  the  heat  conduction  equation 
ar/at  =  (/</pCm)V2r  comes  that  the  rate  of  change  of  energy  is 


d E 
dt 


J  pcJ^AV  =  J  kV2TdV  =  J  kn  -VTdS 

B  B  dB 

-  J  kP(T  -  T„  )dS  =  -kPA(T  -  T„ ) 

dB 


(33) 


where  T  is  the  average  temperature  over  the  surface  66,  that  is 


(34) 


^(x)  #  -|6csin /<(?  -  1)  +  /<cos/<(?  -  1)  for  k^0 


Noting  the  similarity  between  Eq.  (33)  and  the  heat  convection 
equation  in  the  form 


and 

=  ftsin  J(f  +  1)  +  /cos/(?  +  1) 

then  for  k  ¥=  1  and  both  k  and  /  satisfying  Eq.  (27)  we  have  that 

L 

J  tiWM*) dx  =  0  (28) 

-L 

and  for  k  =  l  ¥=  0  we  have  that 


Qconv  —  ~  —hA(T  —  T oo )  (35) 

where  h  is  the  heat  convection  coefficient  [W/m2  I<],  it  is  clear  that 


(36) 


The  boundary  condition  needs  to  assure  continuity  of  temper¬ 
ature  at  the  two  sides  of  the  TEGs.  Therefore,  depending  on  if  the 
two  sides  are  in  thermal  contact  with  a  solid  or  with  a  liquid/gas,  it 
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Fig.  3.  The  system  used  for  all  the  experiments  and  simulations  (not  in  scale). 


will  be  necessary  to  substitute  in  Eq.  (36)  the  appropriate  values, 
where  for  a  solid  h  is 


u  .  _  ,vsohd 

sohd  -  thickness 


(37) 


Fig.  3  illustrates  the  system  used  for  all  the  experiments  and 
simulations  in  this  paper.  It  is  composed  by  a  TEG  sandwiched 
between  a  heater  and  a  water-cooling  block  through  two 
Aluminium  heat-spreading  plates  of  thickness  Lc  =  1  cm  and 
Lh  =  4  cm,  therefore  the  h  coefficients  of  Eq.  (37)  are: 


237  W/mK 
0.01  m 


23,700  W/m2I< 


hH  =  5925  W/m2I< 

The  TEG  is  from  Custom  Thermoelectric  (code:  1261G-7L31- 
04CQ)  and  its  parameters  are: 

Area  :  A  =  0.042m2 
Thickness  :  LTEg  =  0.004  m 
Thermal  Diffusivity :  £~10_7m2/s 
Internal  Resistance  :  Rint~2.5  Q 


It  is  preferable  to  upgrade  Eq.  (36)  to 

_  ^medium 
^TEG 


1 

m 


(38) 


3.  Experimental 

The  solution  worked  out  in  the  previous  section  is  finally  pro¬ 
grammed  in  Matlab,  because  of  the  need  of  a  numerical  calculator 
to  solve  the  Fourier  series  of  Eq.  (21).  Matlab’s  symbolic  toolbox  has 
been  used. 

First,  an  interval  bisection  algorithm  is  used  to  find  the  roots  of 
Eq.  (24),  setting  a  limit  to  the  maximum  value  of  k.  Then  for  every  k 
the  correspondent  coefficient  A'n  is  computed  as  from  Eq.  (30)  and 
finally  the  solution  for  6  is  computed  as  in  Eq.  (31 )  and  scaled  back 
to  temperature.  However  this  solution  does  not  take  into  account 
the  electrical  behaviour  of  the  TE  device  (Seebeck  effect)  and  the 
Peltier  effect,  which  will  therefore  need  to  be  taken  into  account  in 
the  Matlab  model,  as  described  in  this  section. 


The  thermal  diffusivity  of  Bi2Te3  is  in  the  order  [13]  of  10  6,  but 
the  real  one  is  difficult  to  estimate  because  of  the  ceramic  layer 
included  in  the  TEG  for  isolation  and  all  the  interfaces,  which  slow 
down  the  thermal  wave,  therefore  a  value  of  10  7  has  been  used. 

The  internal  resistance  is  a  function  of  the  average  temperature 
of  the  module,  while  the  output  current  and  voltage  depend  on  the 
temperature  gradient.  Following  a  similar  approach  to  the  one 
described  in  a  paper  by  Woo  et  al.  [14],  the  real  curves  of  the 
electrical  characteristic  of  a  TEG  can  be  fitted;  in  this  way,  if  the 
temperature  gradient  AT  and  the  load  current  Iout  are  known,  it  is 
possible  to  mathematically  calculate  the  internal  resistance  R/nt,  the 
output  voltage  Vout  and  the  Seebeck  coefficient  a.  For  the  TEG  used: 

Rint  =  0.00531  AT  +  1.572 

Vout  =  -R-nt/out  +  0.0453AT-0.16 

a  =  0.042  +  0.00263  ^1  -  e~ ^ 

Fig.  4  displays  both  the  experimental  and  the  mathematical  data 
for  V0ut  versus  /out  plotted  for  different  values  of  AT.  As  it  can  be 
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Fig.  4.  Experimental  and  mathematical  data  for  the  electrical  characterisation  of  the  TEG  used,  plotted  for  different  values  of  AT. 
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seen  the  mathematical  functions  well  approximate  the  real 
behaviour  of  the  TEG,  hence  they  can  be  used  in  the  simulation. 

The  thermal  conduction  coefficient  of  the  TEG  /<teg  has  been 
calculated  letting  the  system  get  to  steady-state  with  constant 
input  power  to  the  hot  side  and  open-circuit  at  the  TEG  terminals: 

kTEG  =  P,nA(JH  -  Tc)  ~  1 .4  W/mK  (42) 

^TEG 

in  which  the  power  through  the  TEG  (taking  into  account  the  losses 
of  the  system),  the  temperatures  at  the  sides,  the  area  and  the 
thickness  are  known. 

The  Peltier  effect  can  be  considered  a  parasitic  effect  in  TE 
energy  generation,  because  as  soon  as  some  current  is  drawn  from 
the  TEG,  the  Peltier  effect  brings  power  from  the  hot  to  the  cold 
side,  thus  increasing  the  overall  thermal  conductivity  and  reducing 
the  temperature  gradient  (if  the  heat  powers  to  the  hot  side  and 
from  the  cold  side  remain  the  same).  The  Peltier  effect  is  usually 
greater  than  the  Joule  heating  effect,  therefore  it  is  important  to 
take  it  into  account  in  the  Matlab  model.  In  order  to  do  this,  the 
power  from  the  heater  to  the  TEG  through  the  Aluminium  block  can 
be  considered  as  the  sum  of  the  Peltier  power  and  the  power  from 
Th  oo  to  Th  in  the  heat  equation.  In  this  way  the  solution  is  computed 
in  a  loop.  At  every  iteration  the  heat  equation  is  calculated  and  the 
real  h  coefficient  of  Eq.  (8)  (called  hHreai)  *s  modified  depending  on 
the  Peltier  power  and  called  /iHmocP 

^Hmod  =  ^Hreal  “  a(jHx Tff)  (43) 

Similar  considerations  can  be  applied  to  the  cold  side: 


Heat  equation  with  constant-temperature  boundaries 


Thickness  (mm) 

Fig.  5.  Heat  equation  with  constant  temperature  boundary  conditions,  for  different 
current  loads.  The  hot  side  is  at  x  =  0  mm  while  the  cold  side  is  at  x  =  4  mm. 
Tc  =  15  °C,Th  =  250  °C.  Regardless  of  the  load  current,  the  side  temperatures  are  fixed. 


*Cmod 


aTcI 

Creal  A(Tc_Tcoo) 


(44) 


Therefore  at  the  end  of  every  iteration  fa  and  fa  are  modified 
accordingly  so  that  the  new  steady-state  and  transient  values  will 
account  for  the  Peltier  term.  The  initial  condition  remains  the  same. 


4.  Results  and  discussion 

The  solution  presented  in  Section  2  has  two  main  benefits:  (1 )  it 
allows  dynamic  exchange  of  heat  power  throughout  the  sides  of  the 
TEG  and  (2)  it  allows  the  simulation  of  transients  in  TE  devices. 
These  features  are  analyzed  in  this  section. 

If  the  temperatures  at  the  sides  were  to  be  assumed  constant  as 
in  Eq.  (3),  the  effect  of  Joule  heating  would  be  visible  only  on  the 
points  inside  the  TEG.  Fig.  5  shows  such  a  case  in  a  simulation  in 
which  the  load  current  is  changed.  As  TH  and  Tc  are  constant,  there 
is  no  difference  in  heat  transfer  through  the  TEG  unless  /<teg  is 
modified  to  account  for  the  change  in  internal  Joule  heating. 

On  the  contrary,  the  steady  state  solution  of  Eq.  (16)  is  different 
because  when  a  load  variation  takes  place  TH  and  Tc  can  change 
even  if  Tho 0  and  Tc  oo  remain  constant.  In  this  way  it  is  possible  to 
understand  the  effect  that  different  loads  have  on  the  overall 
conduction  coefficient,  i.e.  to  understand  how  the  temperature 
gradient  across  the  TEG  will  change.  This  can  be  seen  in  Fig.  6,  and  it 
means  that  the  overall  thermal  conductivity  of  the  TEG  changes. 

The  effect  is  more  relevant  at  the  hot  side  because  Th  is  more 
distant  from  THoo  than  Tc  is  from  Tc°°-  From  0  A  to  4  A,  TH  decreases 


Heat  Equation  with  Newton  boundaries 


Fig.  6.  Heat  equation  with  Newton  boundary  conditions  and  Peltier  effect,  for  different 
loads.  The  hot  side  is  at  x  =  -2  mm  while  the  cold  side  is  at  x  =  2  mm. 
TCoo  =  15  °C;Thoo  =  250  °C.  The  hot  side  temperatures  vary  with  the  load  current  (points: 
196,  200,  204.4,  215);  the  cold  side  temperatures  vary  as  well  (points:  18, 19,  20,  24). 
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because  the  Peltier  effect  is  more  relevant  than  the  Joule  heating, 
but  when  the  current  is  8  A  the  Joule  heating  makes  TH  increase 
even  more  than  Thoo  ;  this  would  suggest  an  operating  point  of  great 
efficiency  for  TEGs,  but  unfortunately  it  is  not  physically  possible  to 
make  one  TEG  produce  as  much  current.  Anyway  it  can  clearly  be 
understood  that  the  Peltier  effect  decreases  TH  and  increases  Tc, 
while  the  Joule  heating  increases  both  TH  and  Tc  and  the  temper¬ 
atures  of  the  other  points  inside  the  TEG;  this  effect  depends  on  the 
square  of  the  load  current,  so  it  becomes  quite  relevant  especially  in 
TE  coolers  or  heaters. 

The  fact  that,  when  the  load  current  is  increased,  the  term  aTHI 
of  Peltier  cooling  at  the  hot  side  removes  more  heat  than  the  one 
brought  by  Joule  heating  sets  a  basilar  difficulty  in  proving  with 
experimental  results  the  heat  equation  on  its  own,  because  there 
is  always  the  Peltier  effect  at  the  interface  of  the  pellets  and  an 
ohmic  contact;  nevertheless,  including  the  Peltier  term  into  the 
Matlab  equation  (Eq.  (43)  and  Eq.  (44))  and  making  it  run  itera¬ 
tively  can  provide  satisfactory  simulations  also  for  real  transient 
experiments. 

Steady-state  operating  points  have  been  compared  to  the 
physical  system  described  in  Section  3.  The  system  has  been 
brought  to  steady-state  with  a  load  current  of  0.1  A,  then  changed 
to  1.23  A,  always  with  the  same  input  power  to  the  heater.  The 
comparison  is  shown  in  Table  1  and  shows  that  the  model  can 
predict  Th  and  Tc  quite  accurately. 

Fig.  7  simulates  the  transient  from  an  initial  steady-state 
condition  without  Joule  heating  (open-circuit  load);  at  time  0  s 
the  load  is  suddenly  changed  to  a  current  of  2.25  A.  The  plotted 
curves  show  the  real  and  simulated  transients  of  TH  and  Tc  for  23  s. 
In  the  experiment,  Th oo  is  kept  at  215  °C  with  a  PID  control  on  the 
heater,  while  TCo 0  is  influenced  by  a  1  °C  hysteresis  on  the  water 
chiller.  During  the  experiment  the  cold  temperature  decreases  after 


Table  1 

Experimental  and  simulated  steady-state  operating  points.  THoo  and  TCoo  of  the 
simulation  are  set  as  those  of  the  experiment. 


Load 

Th* o 

Th 

Tc 

TCoo 

Exp  @  0.1  A 

134  °C 

125.4  °C 

22.9  °C 

21  °C 

Exp  @1.23  A 

123  °C 

113  °C 

23.7  °C 

22  °C 

Sim  @  0.1  A 

134  °C 

126.3  °C 

22.9  °C 

21  °C 

Sim  @  1.23  A 

123  °C 

113.8  °C 

23.2  °C 

22  °C 

a  few  seconds  because  the  chiller  starts  cooling  down  to  oppose  the 
change  in  the  water’s  temperature. 

It  can  be  seen  that  the  simulation  can  predict  quite  accurately 
the  evolution  of  Th  and  Tc  during  the  transient,  with  an  error  or 
around  0.5  °C  on  TH. 

The  solution  described  by  Eq.  (31)  could  be  included  into 
a  model  which  takes  into  account  also  the  thermal  masses  of  the 
system,  thus  enabling  the  possibility  of  simulating  both  transients 
and  the  final  operating  point  of  a  complex  system.  This  would 
constitute  a  great  benefit  to  the  actual  engineering  task,  because 
when  designing  a  TE  system  it  is  difficult  to  understand  what  the 
final  operating  point  will  be,  because  the  effects  taking  place  into 
the  TE  device  will  make  the  temperatures  at  the  side  change.  In  fact 
Th oo  and  7c oo  cannot  usually  be  considered  constant.  In  the  case  of 
TE  energy  generation,  the  insertion  of  a  TEG  changes  massively  the 
initially-available  hot  and  cold  temperatures,  depending  on  the 
load  and  on  the  capacity  of  the  system  to  provide  and  remove  hot 
and  cold  power.  As  the  solution  proposed  in  this  paper  analytically 
dynamically  describes  what  happens  inside  the  TE  device,  it 
constitutes  the  starting  point  for  creating  an  accurate  simulation 
model. 


Load  Transient  Experiment  and  Simulation 
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Fig.  7.  Load  change  transient  experiment  and  simulation  results.  /ioad  =  0  A  @  time  =  0  s;  /ioad  =  2.25  A  @  time  >0  s. 
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5.  Conclusions 

This  paper  provides  the  solution  to  the  heat  conduction  equa¬ 
tion  for  TE  devices,  with  internal  Joule  heat  generation,  and 
dynamic  exchange  of  heat  through  the  hot  and  cold  sides.  It  can  be 
used  to  study  the  transient  behaviour  of  TE  devices  as  well  as 
steady-state  operation. 

The  insertion  of  this  solution  into  a  model  which  describes  also 
the  Peltier  effect  allows  the  comparison  of  the  simulated  results  to 
experimental  data,  confirming  the  correctness  of  the  solution. 

Future  work  is  adding  all  the  heat  masses  of  the  system  and  the 
heat  and  electrical  power  flow.  In  this  way  it  would  be  possible  to 
precisely  study  transients  in  TE  systems  and  to  work  out  the  final 
operating  point  when  designing  a  new  system. 
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